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A NEW MATRIX THEOREM AND ITS APPLICATION FOR 
ESTABLISHING INDEPENDENT COORDINATES FOR COMPLEX 
DYNAMICAL SYSTEMS WITH CONSTRAINTS 

By William C. Walton, Jr., and Earl C. Steeves 
Langley Research Center 

SUMMARY 

A new method is presented by which equations of motion of a linear mechanical 
system can be derived in terms of independent coordinates when the system is described 
in terms of coordinates which are not independent but instead are governed by linear 
homogeneous equations of constraint. There is a discussion of the origin in practical 
vibrations analysis of dynamical systems involving equations of constraint. Methods 
previously used for handling such systems are discussed and the new method is demon- 
strated to have the following advantages: (1) For the most general constraint equations, 
solution of the equations is reduced in substance to computing the eigenvalues and eigen- 
vectors of a symmetric matrix; and (2) the method is applicable when there are redun- 
dancies in the equations of constraint. 

INTRODUCTION 

The purpose of this paper is to present a method by which equations of motion of a 
linear mechanical system can be derived in terms of independent coordinates when basic 
information about the system is available in terms of coordinates which are not inde- 
pendent but instead are governed by linear homogeneous equations of constraint. Neces- 
sity for this derivation occurs frequently in practical vibration analysis. It arises 
naturally in studies of the motions of bodies composed of components which have been 
idealized as separate bodies. Experience in analyses of vibrations of engineering struc- 
tures has convinced the authors that this method often offers decided advantages in prac- 
tical computation over methods previously used. 

The method is based on a mathematical theorem designated the "zero eigenvalues 
theorem" which allows the computational procedures to be systematically developed. A 
search of the mathematical literature has been made and nowhere has this result been 
found. 


The paper begins with a background note on dynamical systems involving constraint 
equations. A brief discussion of approaches previously taken in treating such systems 
follows. The zero eigenvalues theorem is proved, and the method of this paper is dis- 
cussed. There is a development of the relationship between the result obtained by the 
method of this paper and the result obtained by the method generally taught in engineering 
textbooks. Two examples of application of the theorem to problems from vibration anal- 
ysis are presented and the numerical considerations involved in practical computing with 
the method are discussed. 


SYMBOLS 

A an R x R partitional submatrix of matrix C 

B an R x (P - R) partitional submatrix of matrix C 

C a constant matrix of order R x P 

D a constant matrix defined by equation (23) 

E a constant matrix defined by equation (14) 

F n elements of vector Q in first example 

G number of positive finite elements of diagonal matrix \ 

H any nonsingular matrix of order P - G 

I identity matrix 

K stiffness matrix referred to coordinates q 

K stiffness matrix referred to coordinates q 

L Lagrangian 

l length of cylinder 

M mass matrix referred to coordinates q 


2 


M mass matrix referred to coordinates q 

m axial wave number 

N number of elements in vector q 

n subscript denoting general element in vector q 

P number of elements in vector q 

p subscript denoting general element in vector q 

Q generalized forces referred to coordinates q 

Q generalized forces referred to coordinates q 

q a vector whose elements are independent coordinates 

q a vector whose elements are dependent coordinates 

q( a ) a partition of q containing R elements 

q(k) a partition of q containing those elements not in q^ a ^ 

q m Q coordinates associated with axisymmetric circumferential 

harmonic 

q coordinate associated with fourth circumferential harmonic 

m,4 

R number of rows in matrix C 

r radius of cylinder 

S quadratic form defined by equation (19) 

T a constant matrix in equation (33c) relating the coordinates q to 

the coordinates q 

T matrix defined by equation (35) 
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U modal matrix of matrix E 

u longitudinal displacement of shell 

v an arbitrary vector of order P 

v a vector defined by equation (22) 

W work of external forces 

x n elements of vector q for first example 

z a vector defined by equation (30) 

z«) a vector whose elements are first G elements of z 

z c displacement of center of disk in £ direction 

«£ rotation about £-axis 

a rotation about 77-axis 

/3 a constant matrix in equation (9) relating dependent coordinates q 

to independent coordinates 

6 variational operator 

4,77,5 coordinates used in second example 

X a real diagonal matrix 

X a real diagonal matrix whose elements are positive elements of X 

Xp pth diagonal element of matrix X 

(') denotes differentiation with respect to time 

( )’ denotes transpose of a matrix 
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denotes a row matrix 


L J 
□ 

O 


denotes a rectangular matrix 
denotes a column matrix 


BACKGROUND 

In conventional analyses of small forced oscillations of mechanical systems, the 
physical system is idealized so that its configuration at any instant is determined by 
specification of a finite number of independent coordinates q^, q 2 , . • q n , . . q^j. 
Then, with approximations allowable because of the assumed smallness of the oscillations, 
the Lagrangian of the system may be expressed as in reference 1 in the form 

L = |qMq -|q'Kq (1) 

where 

(1) q is a column matrix the elements of which are the coordinates q n 

(2) M and K are constant symmetric matrices of order N 

(3) A prime denotes the transpose of a matrix 

(4) A dot denotes differentiation with respect to time. 

When the Lagrangian has the form shown by equation (1) and the coordinates q n 
are independent, Lagrange's equations of motion of the system have the form (see ref. 1): 

Mq + Kq = Q (2) 

In equation (2) Q is a column matrix with N elements. The elements of Q are 
usually called generalized forces. The generalized forces are determined by the fol- 
lowing requirements: Let 5q be an arbitrary infinitesimal variation of the coordinates 
composing the matrix q. Then the work W done by the forces applied to the system 
when these forces act through the displacements produced by the variation shall be given 
by the equation 

W = Q' 6q (3) 

The generalized forces may be functions of the coordinates q n and/or the time 
explicitly. 

Once the equations of motion are known in the form indicated by equation (2) , there 
is a well-established and very effective body of mathematical theory and computational 
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technique for determining the behavior of the system. Often, however, it is much eas- 
ier to express L and W in terms of coordinates which are not independent but which 
are governed by linear homogeneous equations of constraint. (See ref. 2.) Let 
q 1} q 2 , . • q p , • . ., qp represent such a set of coordinates. The constraint equa- 
tions then take the form 

Cq = 0 (4) 

where q is a column matrix the elements of which are the coordinates qp, and where 
C is a constant matrix which has P columns and is, in general, rectangular. 

In terms of the dependent coordinates qp, the Lagrangian will take the form 

L = | q Mq - | q Kq (5) 

where M and K are symmetric matrices of order P. The work W can be found in 
the form 

W = Q 6q (6) 

where 6q is an arbitrary variation of q compatible with the equations of constraint 
(eq. (4)) and Q is a column matrix with P elements which are functions of the coor- 
dinates qp and/or the time explicitly. 

It is useful to know a systematic procedure by which equations of motion in terms 
of independent coordinates, as in equation (2), can be derived by starting with the 
Lagrangian L and the work W in terms of coordinates governed by homogeneous 
equations of constraint as in equations (5) and (6). The object of this paper is to set forth 
such a procedure, but before doing so, it is appropriate to discuss briefly how the prob- 
lem has been solved previously. 


PREVIOUS METHODS 

In the past, the equations of motion in terms of independent coordinates have been 
determined in two ways: 

(1) Through consideration of particular physical or geometrical aspects of a prob- 
lem, the dependent coordinates qp are chosen to impart a very simple form to the equa- 
tions of constraint, which renders easy and obvious determinations of the independent 
coordinates. 

(2) By using one of many variants of Gauss's classical elimination algorithm, the 
equations of constraint are solved as simultaneous equations; these solutions lead to the 
selection of certain of the coordinates as independent coordinates and the expression of 
the remaining coordinates in terms of those which have been selected to be independent. 
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Under the first category of approaches come, for example, those finite-element 
methods of structural analysis in which the coordinates of a free-body element are dis- 
placements and rotations at juncture points among structural elements. In such anal- 
yses the equations of constraint are equalities among appropriate displacements and 
rotations at junctures and equations in which appropriate displacements and rotations 
are set equal to zero at junctures where there are supposed to be rigid constraints. A 
set of independent coordinates is determined by the simple expedient of using a single 
symbol for each set of displacements and rotations which are equated. (See, for example, 
ref. 3.) This idea is the basis of the now widely used procedure of superimposing stiff- 
ness matrices or mass matrices of structural elements to determine a stiffness matrix 
or a mass matrix of an entire structure composed of the connected elements. 

In order to illustrate some advantages of the method to be presented, the method 
generally taught in engineering textbooks is discussed formally. (See, for example, 
ref. 4.) This method belongs in the second category of approaches. It is assumed (usu- 
ally tacitly) that the rank R of the matrix C is equal to the number of rows in C and 
that, therefore, equation (4) may be written as 

Aq (a) + Bq( b) = 0 (7) 

where 

(1) A is an R X R nonsingular constant matrix the columns of which are R 
distinct columns of C 

(2) B is an R x.(P - R) constant matrix the columns of which are those columns 
of C not included in A 

(3) q^ and q^ are column matrices the elements of which are elements of q 
corresponding to the columns in A and B, respectively. The goal is to establish the 
coordinates in q^ as independent coordinates. 

By renumbering the coordinates <jp, it can be arranged that the first R columns 

of C constitute the matrix A and the last P - R columns of C constitute the 

_/ a \ 

matrix B. Correspondingly, the elements of q v ' would be the first R elements of 
q, and the elements of q^ b \ the last P - R elements of q. For convenience in the 
ensuing discussion, it is assumed that such a rearrangement has been made. However, 
as a practical matter, it is very important to note that in order to actually make a suit- 
able rearrangement, one must be able to identify R linearly independent columns of C. 
This identification may not be easy. 

Since A is nonsingular, an inverse of A exists and is unique. Equation (7) is 
satisfied therefore if, and only if, 

q^ = -A _1 Bq^ (8) 
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( 9 ) 


where A - * is the inverse of A. 
are satisfied if, and only if, 


It follows that the equations of constraint (eqs. (4)) 




where 



B 


0 = 


I 


( 10 ) 


In equation (10) I is an identity matrix of order P - R. Thus, the matrix /3 is a 
P x (P - R) matrix. 

Substitution of equation (9) into equations (5) and (6) gives an expression for the 
Lagrangian L and the work W in terms of independent coordinates and in the forms 
shown by equations (1) and (3), respectively. The components of the expressions are 


q = q( b ) 

(Ha) 

M = /3'M(3 

(lib) 

K = /3’K/3 

(11c) 

Q = /3’Q 

(12) 


where it is to be considered that the substitutions from equation (9) have made the ele- 
ments of the matrix Q functions of the coordinates q n and/or the time explicitly. 

It is noted that the matrices M and K thus derived are symmetric. For empha- 
sis, it is pointed out once more that applicability of this method is restricted to the case 
where the rank R of the matrix C is equal to the number of rows of C and that as a 
practical matter in the application, one is required to identify R linearly independent 
columns of the matrix C. 

It is of interest to consider at this point the situation where contrary to the assump- 
tion made in the foregoing discussion, the rank R of C is less than the number of 
rows of C. It is natural for an analyst in idealizing a physical system to try to specify 
only the minimum number of equations necessary to define the system. In this event the 
rows of C are linearly independent and consequently the rank of C equals the number 
of rows. However, in stress and vibration analyses of engineering structures, experience 
is showing that it is possible to specify inadvertently equations which repeat the content 
of equations or combinations of equations previously written. In fact, it can be of great 
convenience to be able to accept redundant equations as may be seen from one of the 
examples given subsequently. Not much has been written on practical methods for 
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solving systems with dependent equations. However, reference 5 provides a good illus- 
tration of how dependent equations of constraint may arise in practice and also a brief 
discussion of an elimination approach used to solve them. 

ZERO EIGENVALUES THEOREM 

The objective is to prove a theorem that is the foundation of the method of this 
paper. Consider the equation 

Cq = 0 (13) 

where C is a matrix with any number of columns and any number of rows. Let P be 
the number of columns and R the number of rows. 

Let a square matrix E of order P be defined by the equation 

E = C’C (14) 

(It may be noted that the determinant of E is the Gramian of the vectors comprising the 
columns of C. (See ref. 6.)) By transposing both sides of equation (14) and using the 
familiar rule for transposing products of matrices, it follows that 

E' = C^C’)' = E (15) 

Therefore, E, being equal to its own transpose, is symmetric. 

It is a well-known property of symmetric matrices that there exist orthogonal 
matrices U of order P satisfying the following equation: 

U'EU = X (16) 

where X is a real diagonal matrix of order P. By customary usage, an orthogonal 
matrix having this property is called a modal matrix of E, and the numbers occupying 
the main diagonal of X are called eigenvalues of E. To say that U is orthogonal 
means that 

U’U = UU' = I (17) 

where the identity matrix I is of order P. 

Let Xp represent the eigenvalue at the intersection of the pth row and the pth 
column of X. If any modal matrix of E is given, it is easy to construct a modal matrix 
of E so that 

H * x 2 g x 3 2 . . . g Xp (18) 

since the positions of the eigenvalues can be reordered simply by reordering the columns 
of the given modal matrix. Henceforth, in this paper when reference is made to a modal 
matrix, it is to be understood that the columns are ordered so that inequality (18) holds. 
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The convention of the preceding paragraph being understood, it is well-known that 
the eigenvalue matrix X associated with a symmetric matrix E is unique; that is, 
any modal matrix of E when substituted for U in equation (16) produces the same 
matrix X. 

Let v be an arbitrary column matrix with P elements, and let a quadratic 
expression S be defined by the equation 

S = v’Ev (19) 

Note that equation (19) can be written as 

S = v'C’Cv (20) 

or alternatively, by using equations (16) and (17), 

S = v’ UU’ EUU'v = v Xv (21) 

where U is any modal matrix of E and where 

v = U'v (22) 

Equation (20) shows that S cannot be negative for any non null v. Equation (21) shows 
that there exist non null forms of v which will make S negative if, and only if, at 
least one of the eigenvalues Xp is negative. If any choice of v is given, then v given 
by v = Uv will satisfy equation (22). (See eq. (17).) Therefore, if one or more of the 
eigenvalues were negative, there would exist forms of v making S negative and this 
condition would be a contradiction. It follows that the eigenvalues Xp are each positive 
or zero. 

Let an R x P matrix D be defined by the equation 

D = CU (23) 

Then from equations (14) and (16), it follows that 

D’D = X (24) 

Let G be the number of the eigenvalues Xp which are positive. Then the last 
P - G eigenvalues are zero. It follows from equation (24) that D has the partitioned 
form indicated by the equation 



In equation (25) the null matrix 0 is an Rx(P-G) matrix and the matrix D is an 
R x G matrix with mutually orthogonal columns so that 

DD = X (26) 

where X is a diagonal matrix the diagonal elements of which are the G positive eigen- 
values of E as indicated by the equation 
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( 27 ) 



X 


G 


By using equation (17), the equations of constraint (eqs. (13)) may be written as 



CUU’q = 0 

(28) 

or alternatively 


Dz = 0 

(29) 

where 


z = U’q 

(30) 


It is clear from equation (25) that equation (29) is satisfied if the first G elements 
of z are zero whatever the last P - G elements of z may be. 

Premultiplying both sides of equation (29) by d' and substitution of equation (25) 
leads to the equation 

Xz( g ) = 0 (31) 

where z^ is a column matrix the elements of which are the first G elements of z. 
Thus, equation (29) cannot be satisfied unless the first G elements of z are zero. 

Solving equation (30) for q gives the unique solution 

q = Uz (32) 


Thus, the following theorem has been proved: 


Theorem : Consider any set of linear homogeneous equations 

Cq = 0 (33a) 

and let the symmetric matrix E be defined by 

E = C'C (33b) 

The most general solution of the equations may be expressed in the form 

q = Tq (33c) 

where T is a matrix whose columns are the columns of any modal matrix 
of E corresponding to eigenvalues of E which have the value zero, and 
where q is an arbitrary column matrix conformable with T. 
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COMPUTATIONAL PROCEDURE 


With the basic theorem from the preceding section, the following procedural outline 
may be set forth. It is assumed that the generalized forces in the column matrix Q 
are functions of time alone, (if Q is a function of the coordinates qp explicitly, addi- 
tional substitutions will be required which depend on the functional form of Q.) 

Given: 

(1) K and M, both constant symmetric matrices of order P 

(2) C, a constant matrix with P columns and any number of rows 

(3) Q, a column matrix with P elements each of which may be a 
function of time. 

Object: 

(1) To compute a matrix T so that: 

(a) The transformation q = Tq relates the dependent coordinates q 
appearing in equations (5) and (6) to a set of independent coordi- 
nates q suitable for use in equation (2) 

(b) The transformation Q = T’Q produces a matrix Q suitable 
for use in equation (2) 

(2) To compute matrices K and M suitable for use in equations (1) 
and (2). 

Procedure: 

(1) Compute E where E = C'C. Then E will be symmetric of order 
P and positive semidefinite 

(2) Compute a modal matrix U and the eigenvalues Xp of the matrix E 
(where p = l,2,3, . . ., P). This operation is standard at modern 
computing installations and, in fact, is one of the most successful appli- 
cations of digital computers 

(3) Identify the columns of U which correspond to zero eigenvalues. This 
step requires attention because in principle one can fairly question the 
possibility of a rigorous distinction between finite eigenvalues and eigen- 
values having the value zero when, as is normal, there is any roundoff 
error in the process by which the eigenvalues are computed. This point 
is discussed in the section "Comments on Numerical Aspects of 
Computation” 
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(4) Assemble a matrix the columns of which are the columns of U cor- 
responding to the eigenvalues having the value zero. This matrix is. 
the required transformation matrix T. Its dimensions are P by 

P - G where G is the number of positive eigenvalues of E 

(5) Compute K and M by the formulas K = T'KT and M = t’mt. 

Then K and M will be symmetric. 

RELATION TO PREVIOUS METHOD 

Equation (33c) gives the most general solution to the equations of constraint 
(eq. (4)). Since the matrix q appearing in equation (33c) is completely arbitrary, the 
solution can just as well be stated in the form 


q = THq = Tq (34) 

where H is any nonsingular square matrix of order P - G and where 

T = TH (35) 

In order that a matrix T may be written as in equation (35), it is both necessary 
and sufficient that the columns of T constitute a set of linearly independent eigenvectors 
of E corresponding to the eigenvalues of E which have the value zero. The eigenvec- 
tors in T will not, in general, be orthonormal nor even orthogonal. The columns of T 
are orthonormal if, and only if, H is an orthogonal matrix and orthogonal if H is a 
diagonal matrix. Proof of these statements will not be made as they amount merely to a 
formal statement of the basic results of that portion of the theory of matrices which deals 
with repeated eigenvalues of a real symmetric matrix. (See, for example, ref. 7.) A 
connection may now be made between the method of this paper and the textbook method as 
given in reference 4. 

By assuming that the column and coordinate rearrangements leading to equation (7) 
have been carried out, one may write 


E = C’C = 


A’ 

B' 


B 



A’ A 

a'b 

~ 


_J """i 

b'a 

B’B 


(36) 
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It follows that 


1 

A’ A I A’B 
1 


-A -1 B 


n 

r 

B T A ! b t b 

1 

1 


I 


U 


( 37 ) 


where in equation (37) the matrix on the right is a P X (P - R) null matrix. It is clear 
from equation (37) that the columns of /3 are linearly independent eigenvectors of E 
corresponding to P - R eigenvalues having the value zero. The RXR submatrix 
A* A is of rank R. The matrix E is consequently of rank R and possesses no more 
than P - R eigenvalues with value zero. 

Thus, the textbook solution, which is a variant of Gauss's classical elimination 
algorithm, is seen to be a solution of the form of equation (34). 

FIRST EXAMPLE 

In the first example, the method of this paper is applied to derive the equations of 
motion of a simple chain of spring-mass elements. The main intent is to illustrate an 
application of the method. However, some points of general interest arise. 

The system consists of five point masses connected by linear massless springs as 
shown in sketch (1). Each of the masses and each of the spring constants are assumed to 
have unit magnitude. The masses may displace only in the horizontal direction and the 
displacement of the nth mass referred to its undeformed position is denoted by x n . A 

positive value of x n is taken to mean displacement 
to the right, and a negative value, displacement to the 
1 2 3 4 5 left. A horizontal external force F n , positive to the 

Sketch (1).- Spring-mass system. ri S ht and a ^tion of time only, acts upon the nth 

mass. 

The five displacements constitute a set of independent coordinates which determine 
the configuration of the system at any instant; and in terms of these coordinates, it is 
easy to write down directly equations of motion of the system in the form 

Mq + Kq = Q (38) 
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where 




1 0 0 0 0 

0 10 0 0 

M = 0 0 10 0 
0 0 0 1 0 

0 0 0 0 1 


(39b) 


(39c) 



Equation (38) has the form of equation (2). Thus, from a practical point of view, 
the method of this paper is not needed for an analysis of the system since the end result 
of the method, equations of motion in terms of indepen- i'VW'v* *W\/\* |WA* 

dent coordinates, is readily obtained by inspection. 

12345678 

However, since the object is to illustrate the method, "" 


let the system be viewed in a different way as illus- 
trated in sketch (2). There the system of sketch (1) 


Sketch (2).- Cut system. 
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is shown figuratively divided into four parts by cuts at the three inner masses to produce 
an eight-mass system. The half circles represent masses of one-half -unit magnitude. 
The displacement of the pth mass of this cut system is denoted by q p . 

It is assumed that three equations of constraint are imposed on the coordinates; 
namely, 

^2 = ^3 (40a) 


q 4 = q 5 

(40b) 

q 6 = q 7 

(40c) 


Thus the coordinates q p are not independent, and from the simple geometric considera- 
tions involved, it is clear that under these equations of constraint, the systems of 
sketch (1) and sketch (2) are the same. In terms of the coordinates q p , a Lagrangian 
of the system may be expressed by an equation like equation (5) with 


1-1000 
-110 0 0 
001-10 
00-110 
0 0 0 0 1 

0 0 0 0 -1 

0 0 0 0 0 

0 0 0 0 0 


0 0 0 

0 0 0 

0 0 0 

0 0 0 

-10 0 
1 0 0 

0 1 -1 

0 -1 1 


M = 


1 

0 

0 

0 

0 

0 

0 

0 


0 

1/2 

0 

0 

0 

0 

0 

0 


0 

0 

1/2 

0 

0 

0 

0 

0 


0 

0 

0 

1/2 

0 

0 

0 

0 


0 0 0 0 

0 0 0 0 

0 0 0 0 

0 0 0 0 

1/20 00 
01/200 
0 01/20 
0 0 0 1 


(41a) 


(41b) 


The work W for the system may be expressed by an equation following the form 
of equation ( 6 ) with 

Q = [Fj, 0, F 2 , 0, F 3 , 0, F 4 , F 5 J (42) 

T 

It is noted that the form indicated by equation (42) for the matrix Q is not unique. The 
following form, for example, will serve equally well: 
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(43) 


q’ = jFp (1/2)F 2 , (1/2)F 2 , (1/3)F 3 , (2/3) Fg, (1/5)F 4 , (4/5) F 4 , F 5 J 


All that is required is that Q when introduced in equation (6) should yield the work done 
during any displacement consistent with the equations of constraint. 


Equations (40), the equations of constraint, may be put in the form of equation (4) 

with 


C = 


0 1 
0 0 
0 0 


-1 0 0 0 0 0 

0 1-10 0 0 
0 0 0 1 -1 0 


(44) 


The first step in the application of the method is to compute the matrix E defined 
by equation (14). This computation yields 


C'C 


00000000 
0 1 -1 00000 

0 -1 1 00000 

0 0 0 1 -1 0 0 0 

0 0 0 -1 1 0 0 0 

000001 -1 0 
00000 -1 10 
00000000 


(45) 


The matrix U which follows is a modal matrix of the matrix E, as may be easily 
verified by substitution of the matrix into equations (16) and (17). 


0 0 

l /\ l 2 0 

-1/1/2 0 

0 l /\ l 2 

0 -1/1/2 

0 0 

0 0 

0 0 


0 1 0 
0 0 1/V2 

0 0 1/1/2 

0 0 0 
0 0 0 
1/V2 0 0 

-l/\|2 0 0 

0 0 0 


0 0 0 

0 00 

000 
1/ V2 0 0 

1/1/2 0 0 

0 1/ V2 0 

0 1/1/2 0 

0 0 1 


(46) 


The matrix X containing the eigenvalues associated with the modal matrix is 
given by 
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inn nil 


X = 


2 

0 

0 

0 

0 

0 

0 

0 


0 

2 

0 

0 

0 

0 

0 

0 


0 

0 

2 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

0 

0 


0 

0 

0 

0 


0 

0 

0 

0 


0 0 
0 0 


0 

0 


0 

0 


(47) 


The first three eigenvalues are finite and the last five have the value zero. Therefore 
the last five columns of U constitute a suitable transformation matrix T. It follows 
that acceptable independent coordinates for describing any configuration of the system 
consistent with the equations of constraint are five coordinates q n related to the dis- 
placements 


qp by the equation 


V 


~ 1 

0 

0 

0 

0 

q 2 


0 

l/& 

0 

0 

0 

q 3 


0 

l/ff 

0 

0 

0 

q 4 


0 

0 

l/Vff 

0 

0 

k 

0 

0 

l/V2 

0 

0 

q 6 


0 

0 

0 

\/n 

0 

q 7 


0 

0 

0 

1/V5 

0 

q 8j 


0 

0 

0 

0 

1 


S" 

q 2 

q 3 

q 4 


q 5 


(48) 




From the equation Q = T'Q it follows also by using either equation (42) or equation (43) 
that generalized forces suitable for use with the coordinates q n are given by 


"Qi 


r F 1 ^ 

Q 2 


W F2 

<^3 

>=< 

CO 

-1^ 

Q 4 





L 5 J 


( 49 ) 
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Completing the steps in the method gives 



fl 

0 

0 

0 

0 ~ 



0 1/2 

0 

0 

0 


M = T’MT = 

0 

0 

1/2 

0 

0 



0 

0 

0 

1/2 0 



0 

0 

0 

0 

1 

— i 



i 

- 1 / V 2 

0 

0 

0 


-iM 

1 

- 1/2 

0 

0 

T'KT = 

0 

- 1/2 

1 

- 1/2 

0 


0 

0 

- 1/2 

1 

- 1 / V2 


0 

0 

0 

- l / V 2 

1 


(50a) 


(50b) 


Equations (49), (50a), and (50b) give all the quantities necessary for writing the equations 
of motion for the spring-mass system in the form of equation (2). By use of equation (48), 
solutions of the equations giving time histories of the coordinates q n can be transformed 
into time histories of the original coordinates qp. If initial conditions consistent with 
the equations of constraint are given in terms of the coordinates qp, the equations 

q = T'q (51) 

and 


q = Tq 


(52) 


may be used to convert them into initial conditions on the coordinates q n . 

It may be noted that the matrices K, M, and Q given in equations (50b), (50a), 
and (49) are not identical to the corresponding matrices in equations (39) which were 
written down directly from simple physical considerations. Either set of matrices forms 
a valid basis for equations of motion of the system of sketch (1). The difference between 
the matrices arises from the fact that the coordinates q n determined by the method of 
this paper are not related to the coordinates q p in the same way as are the displace- 
ment coordinates x n . Equation (48) shows the relationship between the coordinates q n 
and qp whereas coordinates x n and qp are related by the equation 
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*1 


" 1 

0 

0 

0 

0 


fx ^ 
X 1 

^2 


0 

1 

0 

0 

0 


x 2 

<*3 


0 

1 

0 

0 

0 

< 

x 3 > 

- 4 

> = 

0 

0 

1 

0 

0 


x 4 

«5 

0 

0 

1 

0 

0 


x 5 

<*6 


0 

0 

0 

1 

0 


q 7 


0 

0 

0 

1 

0 


*8 


0 

0 

0 

0 

1 



Equation (53) may be written 


q = Tx 


where 


( 53 ) 


(54) 


T = TH (55) 

in which T is the transformation matrix in equation (48) , derived by the method of this 
paper, and 


H = 


1 

0 

0 

0 

0 


0 0 0 

\J2 0 0 

0 € 0 

0 0 \f2 

0 0 0 


0 

0 

0 

0 

1 


(56) 


Thus the coordinates x n which represent displacements of masses are in the category 
of coordinates discussed in connection with equation (34). The foregoing discussion 
illustrates a feature of the method of this paper which should be recognized by anyone 
using the method; that is, the coordinates q n produced by the method are generally 
abstract in character and do not lend themselves to simple physical interpretations. 


It is instructive to reexamine the matrix C in equation (44) and to think about the 
decisions involved in applying the textbook method. Consider the three pairs of displace- 
ments ^ 2 ^ 3 )’ (^ 4 >^ 5 )’ and (^ 6 ^ 7 ) straddling the cuts in sketch (2). Let triplets of dis- 
placements be formed by taking one and only one displacement from each pair, for 
example, the displacements in any such triplet are taken to make up the 

elements of the column q^ in equation (7), the matrix A formed from the corre- 
sponding columns will be nonsingular and the textbook method will succeed. If the ele- 
ments of q^ are chosen from among the eight coordinates q in any other way, the 
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matrix A will be singular. In applying the textbook method to this simple problem, 
recognition of the combinations of coordinates suitable to form must come about 

either from physical insight or from understanding of linear dependence among the col- 
umns of C. In applying the method of this paper, it is not necessary to think directly 
about the physics involved or about the linear dependence. Instead, the problem becomes 
one of finding a modal matrix of E and identifying the columns associated with eigen- 
values having the value zero. Because of the block-diagonal form of E in this case, it 
was possible by inspection to put down exactly a modal matrix and the eigenvalues of E. 
Therefore, all decisions in the application of the method of this paper could be made eas- 
ily on a purely mathematical basis. 

SECOND EXAMPLE 

The purpose of this second example is to show a condition in which redundancies 
in the equations of constraint arise in a natural way. The mechanical system is shown 
in sketch (3). A cylindrical elastic shell is fixed at one end to an immovable base. At 
the other end, a thin massive rigid disk is attached to the wall of the shell by four pins 
placed at 90° intervals around the circumference. Points in the shell wall are assumed 
to displace only longitudinally. 



By adopting an approximation common in practical vibration analysis, the longi- 
tudinal displacement u of a general point in the shell wall is expressed as a linear com- 
bination of a finite number of displacement functions. The expansion assumed is 
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(57) 


u= I( i 


^111,0 + ^m,4 


cos 


40)sin m ^ 
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where m takes on positive integral values and the summation sign indicates summation 
of the terms corresponding to some finite number of selected values of m. The coeffi- 
cients q and q 4 are functions of time alone and serve as coordinates which 

111 jU 

describe the instantaneous configuration of the shell. 

By assuming small displacements, the instantaneous position of the disk is deter- 
mined by specification of three coordinates z c , a^, and defined as follows: 

(1) z c is the displacement of the center of the disk parallel to the longitudinal 
axis of the shell 


(2) a ^ and otq are small rotations about axis | and r), respectively, as shown 
in sketch (3). 

Equating the displacements of the disk to the displacements of the shell at each of 
the four pins gives 


Z C - - Z(Vo + V4) sin f- 

z c + = £(q mj0 + q m>4 ) sin ! f- 


Z c + M 5=2(Vo + V,4) 3m! f I 

zc - ra„ = ^(q m _ 0 + V ;4 )sin 


(58) 


where r is the radius of the cylinder. If, in the summation on the right, only the terms 
corresponding to m = 1 are retained, the equations may be put in the form 


where 


C = 


Cq = 0 


1 

0 

-1 

-1 

-1 

1 

1 

0 

-1 

-1 

1 

0 

1 

-1 

-1 

1 

-1 

0 

-1 

-1 


(59) 


(60) 
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and 



TCL 


V 




TOly 


*1,0 

^1 4 

V- - 


( 61 ) 


It will be clear on inspection that an attempt to determine independent coordinates 
for this system by a straightforward application of the textbook method must fail because 
any choice of the matrix A will lead to a matrix which has at least two proportional 
columns and which is therefore singular. This difficulty stems from the fact that the 
system of equations is redundant; the redundancy may be demonstrated by adding rows 1 
and 3 of matrix C and subtracting row 2 from the result to produce row 4. 

One way to determine independent coordinates would be to discard the fourth equa- 
tion from the system and apply the textbook method to the first three equations. However, 
this approach requires, in general, the following: 

(1) Recognition in the first place that the system is redundant 

(2) Identification of dependent equations 

(3) Identification of a nonsingular submatrix A after redundant equations are 
discarded. 


For the example problem under consideration, the required understanding of the 
structure of the equations may be gained by inspection. In practical work, however, there 
may be many equations of constraint involving many unknowns , and the coefficients making 
up the matrix C will usually not be small integers. Generally, in such situations, little 
of use can be deduced about the system merely by inspection of the matrix of coefficients. 
Also, one cannot always rely on physical insight to detect and understand redundancies. 
Furthermore, there are considerable theoretical and practical difficulties in making com- 
putational tests for redundancy when there is error, such as roundoff error, in the process 
by which the coefficients of the equations of constraint are generated. (See ref. 5.) 

Proceeding now to apply the method of this paper yields the matrix E as 


E = C’C = 


4 

0 

0 

-4 

-4 


0 0 
2 0 
0 2 
0 0 
0 0 


-4 

0 

0 

4 

4 


-4 

0 

0 

4 

4 


(62) 
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The eigenvalues of E are 


X 1 = 12 ^ 
X 2 = 2 
X 3 = 2 > 

X 4 = 0 i 


( 63 ) 


It may be easily verified that the two columns of the matrix T which follow are ortho- 
normal eigenvectors of E corresponding to the two eigenvalues X 4 and X 5 which 
have the value zero. 


iM 

iM 

0 

0 

0 

0 

-1 M 

1/1/2 

2/\/6 

0 


(64) 


Therefore the system may be described by two independent coordinates and qg 
related to the coordinates in q by the equation 

q = Tq (65) 

As can be seen, direct concern with the number and nature of redundancies in the 
equations of constraint is unnecessary when the method of this paper is used. The prob- 
lem reduces in substance to that of determining a modal matrix of E and identifying 
the columns which correspond to eigenvalues with the value of zero. 


COMMENTS ON NUMERICAL ASPECTS OF COMPUTATION 


In the examples it was possible to, put down exactly the matrix C, to carry out 
exactly the multiplication C'C to produce the matrix E, and to determine exactly the 
eigenvalues of E and orthonormal eigenvectors corresponding to the eigenvalues with 
value zero. In practical work, however, numerical error due to roundoff and/or trunca- 
tion may be introduced at any of these three stages of calculation. The extreme effect 
of such errors would, of course, be complete loss of numerical significance in the digits 
representing the eigenvalues of E and the elements of the eigenvectors of E. In the 
event the computation is subject to serious loss of significance, the matrix C is said 
to be "ill-conditioned" with respect to the computing process used. The best indication 
of ill -conditioning is sensitivity of final results to small changes in the elements of C. 
The authors have applied the method of this paper a number of times in practical vibration 
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analysis and have not encountered a situation in which the matrix C is ill-conditioned. 
From general experience, however, the possibility of ill-conditioning must be anticipated 
whenever simultaneous equations are solved numerically, and the method of this paper 
presents no exception to this statement. When an ill-conditioned system arises, the 
recourse most often open is to increase the number of digits carried in the computation. 

If this procedure is attempted in connection with the method of this paper, it should be 
recognized that it may be necessary to increase the carried significant figures in the 
stage of the calculation in which the elements of C are generated as well as in the 
implementation of the multiplication C’C and in the calculation of the eigenvalues and 
eigenvectors of E. 

Another consequence of numerical error is that finite numbers may be generated 
for eigenvalues of E which would be precisely zero if there were no error in the com- 
puting process. Thus, the question is raised, in principle at least, of the possibility of 
rigorous distinction between finite numbers representing finite eigenvalues of E and 
finite numbers representing eigenvalues of E which are, in fact, zero. In the authors' 
experience this possibility has not proved to be a problem in practice. The authors use 
the threshold Jacobi method (ref. 8) to compute the eigenvalues and a modal matrix of E. 
Approximately 15 significant figures are carried throughout the calculation. With this 
procedure, inspection of the eigenvalues computed for E has always revealed two 
clearly distinguishable sets of numbers, the numbers in one set being many orders of 
magnitude smaller than the numbers in the other set. The set of numbers with relatively 
large magnitudes are regarded as finite eigenvalues, and the remaining numbers are con- 
sidered to be eigenvalues with value zero. 

It is not difficult to show that the number of finite eigenvalues of E is equal to R, 
the rank of C. Frequently, R is known from physical or geometric considerations. 

In particular, it often occurs that one knows that the equations of constraint are linearly 
independent in which case the rank R of C is equal to the number of rows of C. 

Such advance information, of course, enhances confidence in the identification of zero 
eigenvalues. 


RESUME 

The purpose of this paper is to present a new method by which equations of motion 
of a linear mechanical system can be derived in terms of independent coordinates when 
the Lagrangian of the system and the generalized forces are expressed with reference 
to coordinates which are not independent but instead are governed by linear homogeneous 
equations of constraint. 
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As background, there are recalled well-known mathematical forms of the 
Lagrangian and the work statement associated with small oscillations of mechanical 
systems. When the coordinates utilized to develop these expressions are independent, 
Lagrange's method may be applied to determine differential equations of motion of the 
system in a form which has been well studied and is subject to powerful methods of 
solution. However, as a matter of convenience, practical analysis frequently starts 
with coordinates which are dependent by virtue of the imposition of linear homogeneous 
equations of constraint. Thus, it is useful to know a relationship by which the dependent 
coordinates can be transformed into a set of independent coordinates. 

Next, there is a discussion of methods previously used for constructing a trans- 
formation from dependent to independent coordinates. In one category of analysis, the 
dependent coordinates are chosen to impart a very simple form to the equations of con- 
straint so that the transformation may be written from inspection. In a second category 
of approaches based on Gaussian elimination, some of the original dependent coordinates 
are selected to be independent coordinates and those of the original coordinates remaining 
are related to those selected to be independent. This procedure has the drawback that it 
requires in effect the identification of a square nonsingular submatrix in the matrix of the 
coefficients of the equations of constraint. 

The third part of the paper is devoted to the basic result, a theorem which serves 
as a foundation for a new method for constructing the transformation from independent to 
dependent coordinates. 

Theorem: Let a real symmetric matrix be constructed by multiplying the 
matrix of coefficients in the equations of constraint by the transpose of the 
same matrix. Consider any modal matrix of the symmetric matrix so defined 
and select from the modal matrix the columns corresponding to eigenvalues 
with value zero. Then the matrix of these columns is a legitimate transfor- 
mation relating the original dependent coordinates to a set of independent 
coordinates. 

The advantages of constructing the transformation matrix by this method are (1) compu- 
tation is reduced essentially to generating a modal matrix and the eigenvalues of a real 
symmetric matrix and (2) the method is applicable to systems where there are redundant 
equations among the equations of constraint. 

Computing procedures for applying the method are given in outline form, the method 
is applied to two simple physical problems, and numerical considerations in the applica- 
tion of the method are discussed. Also a connection is made between the method of this 
paper and the method generally given in engineering textbooks. In addition to illustrating 
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the method, the examples bring out the abstract nature of the coordinates produced by 
the method and indicate how redundancies in the equations of constraint may arise in 
practical vibration analyses. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va., August 11, 1969. 
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